import random
import numpy as np
import math
from matplotlib import pyplot
import gamma

def loadSizeDistribution():
    name = '/home/jonk/Documentos/papers/Bombard/HQIronParticleSizeDistribution.dat'
    size = [0.0]
    freq = [0.0]
    for line in file(name):
        line = line.replace('\n','')
        val=line.split(' ')
        res=[]
        for v in val:
            if v!='':
                res.append(float(v))
        size.append(res[0])
        freq.append(res[1])
    #Normalizacion
    suma=0.0
    for i in range(1,len(freq)):
        suma += (size[i]-size[i-1])*(freq[i]+freq[i-1])*0.5
    freqnorm=[]
    for f in freq:
        freqnorm.append(f/suma)
    return [size,freqnorm]

def generateRandom(alpha,xm,npoint):
    kappa=alpha+1
    theta=xm/(alpha+1)
    values=[]
    for i in range(npoint):
        values.append(random.gammavariate(kappa,theta))
    return values;

def histogram(values, nhist):
    maxValue=max(values)
    minValue=min(values)
    suma=0.0
    
    dist = np.zeros((nhist+1))
    count = np.zeros((nhist+1))
    step=(maxValue-minValue)/float(nhist)
    for val in values:
        dnd=int(math.floor((val-minValue)/step))
        dist[dnd]+=val
        count[dnd]+=1.0
    
    #Normalizacion
    for i in range(1,len(count)):
        suma+=step*(count[i]+count[i-1])*0.5
    x=[]
    y=[]
    for i in range(len(count)):
        x.append(minValue+float(i)*step)
        y.append(count[i]/suma)
    return [x,y]

def schultz(x,xm,alpha):
    C1=(alpha+1.0)/xm
    uno = math.pow(C1,alpha+1)/gamma.Gamma(alpha+1)
    fun = uno*math.pow(x,alpha)*math.exp(-C1*x)
    return fun

def testRandomSchultz(alpha,xm,nrandom):
    values = generateRandom(alpha,xm,npoint)
    nhist=100
    [x,y] = histogram(values, nhist)
    pyplot.plot(x,y,'ro')
    desde=min(x)
    hasta=max(x)
    puntos=nhist
    paso =(hasta-desde)/float(puntos)
    x=[]
    y=[]
    for i in range(puntos):
        xaux=(float(i))*paso
        x.append(xaux)
        y.append(schultz(xaux,xm,alpha))
    pyplot.plot(x,y,'b-')
    pyplot.show()

def testRandomIronHQ():
    xm=0.8
    alpha =0.4
    nrandom=1000000
    values = generateRandom(alpha,xm,npoint)
    values2=[]
    for v in values:
        values2.append(v+0.3)
    vmedio = np.mean(values2)
    print "Diametro medio:",vmedio
    nhist=1000
    [x,y] = histogram(values2, nhist)
    pyplot.plot(x,y,'ro')
    [size,freq] = loadSizeDistribution()
    pyplot.plot(size,freq,'b-')
    pyplot.show()

        
xm=0.8
alpha =0.4
npoint=1000
#testRandomIronHQ()
values = generateRandom(alpha,xm,npoint)
values2=[]
for v in values:
    values2.append(v/1.1)
f = open("/home/jonk/projects/GenerateConfig/diametros.dat",'w')
for v in values2:
    f.write(str(v)+'\n')
f.close()